# -*- coding: utf-8 -*-
"""
波士顿房价预测 - 线性回归模型
结合数据探索、模型评估和可视化分析
"""
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.metrics import explained_variance_score, median_absolute_error
from sklearn.preprocessing import StandardScaler

# 设置中文字体和图形样式
plt.rcParams['font.sans-serif'] = ['SimHei']  # 支持中文显示
plt.rcParams['axes.unicode_minus'] = False    # 正确显示负号
column_names = [
    'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE',
    'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV'
]
def load_boston_data():
    """
    从本地CSV文件加载波士顿房价数据集
    返回: 特征数据框和目标序列
    """
    # 确保data目录存在（由main函数创建）
    file_path = './data/boston_housing.csv'
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"数据文件 {file_path} 不存在，请检查目录结构")
    
    df = pd.read_csv(file_path, names=column_names)
    if 'MEDV' not in df.columns:
        raise ValueError("数据集必须包含'MEDV'列")
    
    # 数据清洗：移除极端异常值（处理明显不合理的房价）
    Q1 = df['MEDV'].quantile(0.01)
    Q3 = df['MEDV'].quantile(0.99)
    df = df[(df['MEDV'] >= Q1) & (df['MEDV'] <= Q3)]
    
    feature_names = df.columns.tolist()
    feature_names.remove('MEDV')
    return df, feature_names
def explore_data(df):
    """
    数据探索分析
    """
    print("=" * 50)
    print("数据探索分析")
    print("=" * 50)
    print(f"数据集维度: {df.shape}")
    print(f"\n前5行数据:")
    print(df.head())
    print(f"\n数据基本信息:")
    print(df.info())
    print(f"\n描述性统计:")
    print(df.describe())
    
    # 检查缺失值
    missing_values = df.isnull().sum()
    if missing_values.any():
        print(f"\n缺失值情况:\n{missing_values[missing_values > 0]}")
    else:
        print("\n无缺失值")
    
    # 特征相关性分析
    print(f"\n目标变量MEDV与各特征的相关性:")
    correlations = df.corr()['MEDV'].sort_values(ascending=False)
    print(correlations)

def visualize_data(df):
    """
    数据可视化分析
    """
    fig, axes = plt.subplots(2, 2, figsize=(15, 12))
    
    # 房价分布直方图
    axes[0, 0].hist(df['MEDV'], bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    axes[0, 0].set_xlabel('房屋价格（千美元）')
    axes[0, 0].set_ylabel('频数')
    axes[0, 0].set_title('房价分布直方图')
    axes[0, 0].grid(True, alpha=0.3)
    

    
    # 关键特征与房价的关系散点图（选择两个重要特征）
    axes[1, 0].scatter(df['RM'], df['MEDV'], alpha=0.6, color='green')
    axes[1, 0].set_xlabel('房间数量 (RM)')
    axes[1, 0].set_ylabel('价格')
    axes[1, 0].set_title('房间数量 vs 房价')
    axes[1, 0].grid(True, alpha=0.3)
    
    axes[1, 1].scatter(df['LSTAT'], df['MEDV'], alpha=0.6, color='red')
    axes[1, 1].set_xlabel('低收入人口比例 (LSTAT)')
    axes[1, 1].set_ylabel('价格')
    axes[1, 1].set_title('低收入人口比例 vs 房价')
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

def build_and_evaluate_model(X, y):
    """
    构建线性回归模型并进行评估
    返回: 训练好的模型和评估结果
    """
    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42
    )
    
    # 数据标准化，提升线性回归性能
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train)
    X_test_scaled = scaler.transform(X_test)
    
    # 构建并训练线性回归模型
    model = LinearRegression()
    model.fit(X_train_scaled, y_train)
    
    # 预测
    y_pred = model.predict(X_test_scaled)
    
    # 计算评估指标
    metrics = {
        'MSE': mean_squared_error(y_test, y_pred),
        'RMSE': np.sqrt(mean_squared_error(y_test, y_pred)),
        'MAE': mean_absolute_error(y_test, y_pred),
        'R2': r2_score(y_test, y_pred),
        '解释方差': explained_variance_score(y_test, y_pred),
        '中值绝对误差': median_absolute_error(y_test, y_pred)
    }
    
    # 交叉验证评估
    cv_scores = cross_val_score(model, X_train_scaled, y_train, cv=5, scoring='r2')
    metrics['交叉验证R2均值'] = cv_scores.mean()
    metrics['交叉验证R2标准差'] = cv_scores.std()
    
    return model, metrics, y_test, y_pred, X_train, X_test, scaler

def print_model_results(model, metrics, feature_names):
    """
    打印模型结果和评估指标
    """
    print("=" * 50)
    print("模型评估结果")
    print("=" * 50)
    
    print(f"截距项: {model.intercept_:.4f}")
    print("\n特征系数:")
    for name, coef in zip(feature_names, model.coef_):
        print(f"  {name}: {coef:>8.4f}")
    
    print("\n评估指标:")
    for metric, value in metrics.items():
        if metric == 'R2':
            print(f"  {metric}: {value:.4f}")
        else:
            print(f"  {metric}: {value:.4f}")

def visualize_results(y_test, y_pred):
    """
    可视化模型预测结果
    """
    fig, axes = plt.subplots(1, 2, figsize=(15, 6))
    
    # 真实值 vs 预测值散点图
    axes[0].scatter(y_test, y_pred, alpha=0.7, color='blue')
    axes[0].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
    axes[0].set_xlabel('真实价格')
    axes[0].set_ylabel('预测价格')
    axes[0].set_title('真实值 vs 预测值')
    axes[0].grid(True, alpha=0.3)
    
    # 添加R²到图中
    r2 = r2_score(y_test, y_pred)
    axes[0].text(0.05, 0.95, f'R² = {r2:.4f}', transform=axes[0].transAxes, 
                bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8))
    
    # 预测值与真实值对比折线图（前50个样本）
    indices = range(min(50, len(y_test)))
    axes[1].plot(indices, y_test.values[:50], 'o-', label='真实值', linewidth=2, markersize=6)
    axes[1].plot(indices, y_pred[:50], 's-', label='预测值', linewidth=2, markersize=4)
    axes[1].set_xlabel('样本索引')
    axes[1].set_ylabel('价格')
    axes[1].set_title('预测值与真实值对比（前50个样本）')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

def main():
    print("波士顿房价预测分析")
    print("=" * 50)
    
    # 1. 加载数据
    df, feature_names = load_boston_data()
    X = df.drop('MEDV', axis=1)
    y = df['MEDV']
    
    # 2. 数据探索
    explore_data(df)
    
    # 3. 数据可视化
    visualize_data(df)
    
    # 4. 构建和评估模型
    model, metrics, y_test, y_pred, X_train, X_test, scaler = build_and_evaluate_model(X, y)
    
    # 5. 打印结果
    print_model_results(model, metrics, feature_names)
    
    # 6. 结果可视化
    visualize_results(y_test, y_pred)
    
    print("\n分析完成！")


if __name__ == "__main__":
    # 创建目录
    os.makedirs('./data', exist_ok=True)
    main()